Introduction

This file contains essential commands from Chapter 7: Exploratory data analysis of the textbook r4ds and corresponding examples and exercises. A command is considered “essential” when you really need to know it and need to know how to use it to succeed in this course.

All ds4psy essentials so far:

Nr.	Topic
0.	Syllabus
1.	Basic R concepts and commands
2.	Visualizing data
3.	Transforming data
4.	Exploring data (EDA)
+.	Datasets

Course coordinates

Taught at the University of Konstanz by Hansjörg Neth (h.neth@uni.kn, SPDS, office D507).
Winter 2018/2019: Mondays, 13:30–15:00, C511.
Links to current course syllabus | ZeUS | Ilias

Preparations

Create an R script (.R) or an R-Markdown file (.Rmd) and load the R packages of the tidyverse. (Hint: Structure your script by inserting spaces, meaningful comments, and sections.)

## Exploring data (EDA) | ds4psy
## 2018 11 28
## ----------------------------

## Preparations: ----------

library(tidyverse)

## 1. Topic: ----------

# etc.

## End of file (eof). ----------

To use R Markdown, create a corresponding file and save it with the .Rmd extension (e.g., by selecting File > New File > R Markdown). For instructions on combining text and code, see Chapter 27: R Markdown of our textbook, or use one of the following templates:

minimal template: rmd_template_s [in .Rmd | .html format]
medium template: rmd_template_m [in .Rmd | .html format]
explicit explanations: Rmarkdown_basics [in .Rmd | .html format]

Hint: Try to knit your .Rmd file immediately after saving it and marvel at the beauty of the resulting .html-file. If this works, keep doing this routinely from now on, putting all your R-code into code chunks, and any text (like headings or conclusions) that describes or explains what you are doing outside of them. From now on, you can share your .html output files, rather than your .Rmd source files when showing off your R and data science skills.

Exploring data

Introduction

In the following, we refine and use what we have learned so far to explore data. Practically, this session combines what we have learned about ggplot2 (in Chapter 3: Data visualization) and about dplyr (in Chapter 5: Data transformation) to explore datasets.

Key concepts

Important concepts in this session include:

missing values (NA) vs. unusual and extreme values (e.g., outliers),
different types of variables (e.g., continuous vs. categorical variables),
relationships between variables (of different types),
trends (i.e., developments over time), and
many different types of plots.

See Chapter 7: Exploratory data analysis (EDA) and the links provided below for additional information.

What is EDA?

According to Grolemund & Wickham (2017, Chapter 7) exploratory data analysis (EDA) primarily is a state of mind. Rather than following a strict set of rules, EDA is an initial phase of familiarising yourself with a new data set by actively engaging in an iterative cycle during which you

Generate questions about data.
Search for answers by visualising, transforming, and modelling data.
Use answers to refine your questions and/or generate new questions.

Philosophically speaking, EDA is a data scientist’s way of doing hermeneutics (see Wikipedia or The Stanford Encyclopedia of Philosophy for definitions) to get a grip on data.

The goal and purpose of EDA is to gain an overview of a dataset. This includes an idea of its dimensions and the number and types of variables contained in the data, but also a more detailed idea of the distribution of variables, their potential relations to each other, and potential problems (e.g., missing values or outliers).

When using the tools provided by the tidyverse, the fastest way to gain insights into a dataset is a series of dplyr calls and ggplot2 graphs. However, creating good graphs is both an art and a craft. The key to creating good graphs requires answering 2 sets of questions:

Knowing the intended type of plot. This includes answering functional questions like
- What is the goal or purpose of this plot?
- What are possible plot types for this purpose?
- Which of these would be the most appropriate plot here?
Knowing the number and type of variables to be plotted. This includes answering data-related questions like
- How many variables are to be plotted and how are they mapped to dimensions and aesthetics?
- Are these variables categorical or continuous?
- Do some variables control or qualify (e.g., group) the values of others?

Even when these questions are answered, creating beautiful and informative graphs with ggplot requires dedicated practice, experience, and trial-and-error experimentation. In addition, an new dataset is rarely in a condition and shape to directly allow plotting. Instead, we typically have to interleave commands for plotting and transforming data to wrangle variables into specific shapes or types. Thus, calls to ggplot2 usually occur in combination with other tidyverse commands (stemming from dplyr, forcats, tidyr, readr, tibble, etc.).

What needs to be done in any specific case depends on the details of the data and your current goals. Also, there is no pre-defined end to an EDA and no clear boundary between the processes of exploration and the confirmation (or falsification) of expectations. Typically, social scientists use EDA to check and understand a dataset, before using statistics to test specific hypotheses.

While any actual EDA is tailored to specific features of the data and current research goals, this session highlights some common themes that occur in most cases. As we will see, the packages dplyr and ggplot provide a powerful set of tools that allow us to puzzle about data by asking and answering questions. Ideally, you will soon experience that actively engaging in EDA really feels like “doing research”, as it requires us to probe and scrutinize data and follow up on the ideas and hypotheses that we may discover along the way.

To illustrate the process and spirit of EDA, we will use a real and fairly complex data set, which was collected to measure the short- and long-term effects of positive psychology interventions (see Dataset: Positive Psychology for details).¹ And while it’s understandable that you’re mostly focusing on the code and commands, let’s try to keep an eye on the meaning of our data and stay curious about the results that we may find.

Principles and practices

In the following, we will explain some principles that endorse and promote the ideal of transparent data analysis and reproducible research. While such practices are indispensable when working in a team of colleagues and the wider scientific community, organizing your workflow in a more consistent fashion is also beneficial for your other projects and your future self.

Setting the stage

Before embarking on any actual analysis, we should pay tribute to 2 principles that seem so simple that it’s easy to overlook their benefits:

Principle 1: Start with a clean slate and explicitly load all data and all packages required.

Although RStudio IDE provides options for saving your workspace and for loading data or packages by clicking buttons or checking boxes, doing so renders the process of your analysis intransparent for anyone not observing all your actions. To make sure that others and yourself can repeat the same sequence of steps tomorrow or next year, it is advisable to always start with a clean slate and explicitly state which data and packages are required for the current analysis (unless there are good reasons to deviate from them):²

Cleaning up

Clean your workspace and define some custom objects (like colors or functions):

# Housekeeping:
rm(list = ls())  # cleans ALL objects in current R environment (without asking for confirmation)!

# Customizations:
seeblau <- rgb(0, 169, 224, names = "seeblau", maxColorValue = 255)  # seeblau.4 of uni.kn color scheme

# source(file = "my_custom_functions.R")

Loading packages and data

Load all required packages and data files (see Dataset: Positive Psychology for details):

# Load packages:
library(tidyverse)
library(knitr)
library(rmarkdown)

# Load csv-data files (from online links):
# See <http://rpository.com/ds4psy/essentials/datasets.html> for details. 

# 1. Participant data: 
posPsy_p_info <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_participants.csv")

# 2. Original DVs in long format:
AHI_CESD <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_AHI_CESD.csv")

# 4. Transformed and corrected version of all data (in wide format): 
posPsy_wide <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_data_wide.csv")

Clarity

Our second principle is just as innocuous, but becomes crucial as your analyses get longer and more complicated:

Principle 2: Structure and comment your analysis.

While mentioning or adhering to this principle may seem superfluous at first, you and your collaborators will benefit from transparent structure and clear comments as things get longer, messier, and tougher (e.g., when analyses involve different tests or are distributed over multiple datasets and scripts). A consistent document structure, meaningful object names, and clear comments are indispensable when sharing your data or scripts with others. Using R Markdown (see the links and templates above) encourages and facilitates structuring a document, provided that you use some consistent way of combining text and code chunks. Good comments can structure a longer sequences of text or code into parts and briefly state the content of each part (what). More importantly, however, comments should note the goals and reasons for your decisions (i.e., explain why, rather than how you did something).

When trying to be clear, even seemingly trivial and tiny things matter. Some examples include:

Your text should contain clear headings, which includes using consistent levels and lengths, as well as brief paragraphs of text that motivate code chunks and summarize intermediate results and conclusions.
Messy code often works, but is just as annoying as bad spelling: yucAnlivewith outit, but it sure makes life a lot easier. So do yourself and others a favor by writing human-readable code (e.g., putting spaces between operators and objects, putting new steps on new lines, and indenting your code according to the structure of its syntax.³ Incidentally, one of the most powerful ways to increase the clarity of code is by using space: While adding horizontal space seperates symbols and allows structuring things into levels and columns, adding vertical space structures code/text into different steps/sentences, chunks/paragraphs, or larger functional units/sections. And whereas printed media typically need to cut and save characters, space comes for free in electronic documents – so make use of this great resource!

Use clear and consistent punctuation in your comments: A heading or comment that refers to subsequent lines should end with a colon (“:”), whereas a comment or conclusion that refers to something before or above itself should end with a period (“.”).

Yes, all this sounds terribly pedantic and nitpicky, but trust me: As soon as you’re trying to read someone else’s code, you will be grateful if they tried to be clear and consistent.

Safety

It is only human to make occasonal errors. To make sure that errors are not too costly, a good habit is to occasionally save intermediate steps. As we have seen, R allows us to create a copy y of some object x as simply as y <- x. This motivates:

Principle 3: Make copies (and copies of copies) of your data.

Making copies (or spanning some similar safety net) is useful for your objects, data files, and your scripts for analyzing them. Although it also is a good idea to always keep a `current master version'' of a data file, it is advisable to occasionally copy your current data and then work on the copy (e.g., before trying out some unfamiliar analysis or command). As creating copies and then working on them is so easy, it can even save you from remembering and repeatedly typing complicated object names (e.g., plotting some variables ofdfrather thanposPsy_p_info`).

Copies of objects

Most importantly, however, working on a copy of data allways allows recovering the last sound version if things go terribly wrong. Here’s an example:

df <- posPsy_wide  # copy data (and simplify name)
dim(df)            # 295 cases x 294 variables
#> [1] 295 294

sum(is.na(df))  # 37440 missing values
#> [1] 37440
df$new <- NA    # create a new column with NA values
df$new          # check the new column
#>  [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [24] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [47] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [70] NA NA NA NA NA NA
#>  [ reached getOption("max.print") -- omitted 220 entries ]
df     <- NA    # create another column with NA values

# Ooops... 
df  # looks like we accidentally killed our data!  
#> [1] NA

# But do not despair:
df <- posPsy_wide  # Here it is again: 
dim(df)            # 295 cases x 294 variables
#> [1] 295 294

Copies of data

During a long and complicated analysis, it is advisable to save an external copy of your current data when having reached an important intermediate step. While there are many different formats in which data files can be stored, a good option for most files is csv (comma-separated-values), which can easily be read by most humans and machines. Importantly, always verify that the external copy preserves the key information contained in your current original:

# Write out current data (in csv-format):
write_csv(df, path = "my_precious_data.csv")

# Re-read external data (into a different object):
df_2 <- read_csv(file = "my_precious_data.csv")

# Verify that original and copy are identical:
all.equal(df, df_2)
#> [1] TRUE

Screening data

Screening data involves checking the dimensions of data, the types of variables, missing or atypical values, and the validity of observations.

Basic properties

What do you want to know immediately after loading a data file?

Principle 4: Know your data (variables and observations).

This principle is a short version of the following: Directly after loading a new file, you typically want to know

the dimensions of your data,
the types of your variables (columns), and
the semantics of your observations (rows).

Here, we first inspect the original data file, which we read in from http://rpository.com/ds4psy/data/posPsy_AHI_CESD.csv above:

df <- AHI_CESD  # copy the data

dim(df)  # 992 cases x 50 variables
#> [1] 992  50
df <- as_tibble(df)  # (in case df isn't a tibble already)

# Get an initial overview: 
df
#> # A tibble: 992 x 50
#>       id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04 ahi05
#>    <int>    <int>        <dbl>        <int> <int> <int> <int> <int> <int>
#>  1     1        0         0               4     2     3     2     3     3
#>  2     1        1        11.8             4     3     3     4     3     3
#>  3     2        0         0               1     3     4     3     4     2
#>  4     2        1         8.02            1     3     4     4     4     3
#>  5     2        2        14.3             1     3     4     4     4     3
#>  6     2        3        32.0             1     3     4     4     4     4
#>  7     2        4        92.2             1     3     3     2     3     3
#>  8     2        5       182.              1     3     3     3     4     2
#>  9     3        0         0               4     3     3     2     4     2
#> 10     3        2        16.4             4     3     3     3     4     4
#> # ... with 982 more rows, and 41 more variables: ahi06 <int>, ahi07 <int>,
#> #   ahi08 <int>, ahi09 <int>, ahi10 <int>, ahi11 <int>, ahi12 <int>,
#> #   ahi13 <int>, ahi14 <int>, ahi15 <int>, ahi16 <int>, ahi17 <int>,
#> #   ahi18 <int>, ahi19 <int>, ahi20 <int>, ahi21 <int>, ahi22 <int>,
#> #   ahi23 <int>, ahi24 <int>, cesd01 <int>, cesd02 <int>, cesd03 <int>,
#> #   cesd04 <int>, cesd05 <int>, cesd06 <int>, cesd07 <int>, cesd08 <int>,
#> #   cesd09 <int>, cesd10 <int>, cesd11 <int>, cesd12 <int>, cesd13 <int>,
#> #   cesd14 <int>, cesd15 <int>, cesd16 <int>, cesd17 <int>, cesd18 <int>,
#> #   cesd19 <int>, cesd20 <int>, ahiTotal <int>, cesdTotal <int>

## Other ways to probe df: 
# names(df)    # Note variable names
# glimpse(df)  # Note types of variables
# summary(df)  # Note range of values

We note that AHI_CESD contains 992 observations (rows) and 50 variables (columns). Most of the variables are integers and – judging from their names – many belong together (in blocks). Their names and ranges suggest that they stem from questionnaires or scales with fixed answer categories (e.g., from 1 to 5 for ahi__, and from 1 to 4 for cesd__). Overall, a row of data is a measurement of one participant (characterised by its id and an intervention value) at one occasion (characterised by the value of occasion and elapsed.days) on two scales (ahi and cesd).

Unusual values

Principle 5: Know and deal with unusual values.

Unusual values include:

- missing values (`NA` or `-77`, etc.),  
- extreme values (outliers), 
- other unusual values (e.g., unexpected and impossible values).

Missing values

What are missing values? In R, missing values are identified by NA. Note that NA is different from NULL: NULL represents the null object (i.e., something is undefined). By contrast, NA means that some value is absent.

Both NA and NULL are yet to be distinguished from NaN values.

## Checking for NA, NULL, and NaN:

# NA: 
is.na(NA)   # TRUE
#> [1] TRUE
is.na(NULL) # 0 (NULL)
#> logical(0)
is.na(NaN)  # TRUE!
#> [1] TRUE

# NULL:
is.null(NULL) # TRUE
#> [1] TRUE
is.null(NA)   # FALSE
#> [1] FALSE
is.null(NaN)  # FALSE
#> [1] FALSE

# NaN:
is.nan(NaN)  # TRUE 
#> [1] TRUE
is.nan(0/0)  # TRUE
#> [1] TRUE
is.nan(1/0)  # FALSE, as it is +Inf
#> [1] FALSE
is.nan(-1/0) # FALSE, as it is -Inf
#> [1] FALSE
is.nan(0/1)  # FALSE (as it is 0)
#> [1] FALSE
is.nan(NA)   # FALSE
#> [1] FALSE
is.nan(NULL) # 0 (logical)
#> logical(0)

# Note different modes:
all.equal(NULL, NA) # Note: NULL vs. logical
#> [1] "target is NULL, current is logical"
all.equal("", NA)   # Note: character vs. logical
#> [1] "Modes: character, logical"              
#> [2] "target is character, current is logical"
all.equal(NA, NaN)  # Note: logical vs. numeric 
#> [1] "Modes: logical, numeric"              
#> [2] "target is logical, current is numeric"

Missing data values require attention and special treatment. They should be identified, counted and/or visualized, and often removed or recoded (e.g., replaced by other values).

Counting NA values:

df <- AHI_CESD  # copy the data

is.na(df)       # asks is.na() for every value in df 
#>           id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04
#>   [1,] FALSE    FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
#>        ahi05 ahi06 ahi07 ahi08 ahi09 ahi10 ahi11 ahi12 ahi13 ahi14 ahi15
#>   [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#>        ahi16 ahi17 ahi18 ahi19 ahi20 ahi21 ahi22 ahi23 ahi24 cesd01 cesd02
#>   [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  FALSE  FALSE
#>        cesd03 cesd04 cesd05 cesd06 cesd07 cesd08 cesd09 cesd10 cesd11
#>   [1,]  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE
#>        cesd12 cesd13 cesd14 cesd15 cesd16 cesd17 cesd18 cesd19 cesd20
#>   [1,]  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE  FALSE
#>        ahiTotal cesdTotal
#>   [1,]    FALSE     FALSE
#>  [ reached getOption("max.print") -- omitted 991 rows ]
sum(is.na(df))  # sum of all instances of TRUE
#> [1] 0

Since df does not seem to include any NA values, we create a tibble tb that does:

# Create a df with NA values:
set.seed(42) # for replicability
nrows <- 6
ncols <- 6
tb <- as_tibble(matrix(sample(c(1:12, -66, -77), nrows * ncols, replace = TRUE), nrows, ncols)) # create some df
tb[tb > 9] <- NA # SET certain values to NA
# tb

Counting and recoding NA values:

# Count NA values:
sum(is.na(tb))  # => 10 NA values
#> [1] 10

# The function `complete.cases(x)` returns a logical vector indicating which cases in tb are complete:
complete.cases(tb)         # test every row/case for completeness
#> [1] FALSE FALSE FALSE FALSE  TRUE FALSE
sum(complete.cases(tb))    # count cases/rows with complete cases
#> [1] 1
which(complete.cases(tb))  # indices of case(s)/row(s) with complete cases
#> [1] 5

tb[complete.cases(tb), ]  # list all complete rows/cases in tb
#> # A tibble: 1 x 6
#>      V1    V2    V3    V4    V5    V6
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1     9     7   -77   -77     7     1
tb[!complete.cases(tb), ] # list all rows/cases with NA values in tb
#> # A tibble: 5 x 6
#>      V1    V2    V3    V4    V5    V6
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1   -66    NA   -77     7     2    NA
#> 2   -77     2     4     8     8    NA
#> 3     5    NA     7   -66     6     6
#> 4    NA    NA   -77     2   -66    NA
#> 5     8    NA     2   -77    NA    NA

# Recode all instances of NA as -99:
tb
#> # A tibble: 6 x 6
#>      V1    V2    V3    V4    V5    V6
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1   -66    NA   -77     7     2    NA
#> 2   -77     2     4     8     8    NA
#> 3     5    NA     7   -66     6     6
#> 4    NA    NA   -77     2   -66    NA
#> 5     9     7   -77   -77     7     1
#> 6     8    NA     2   -77    NA    NA
tb[is.na(tb)] <- -99  # recode NA values as - 99
tb
#> # A tibble: 6 x 6
#>      V1    V2    V3    V4    V5    V6
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1   -66   -99   -77     7     2   -99
#> 2   -77     2     4     8     8   -99
#> 3     5   -99     7   -66     6     6
#> 4   -99   -99   -77     2   -66   -99
#> 5     9     7   -77   -77     7     1
#> 6     8   -99     2   -77   -99   -99

More frequently, special values (like -66, -77 and -99) indicate missing values. To deal with them in R, we recode all instances of these values in tb as NA:

tb
#> # A tibble: 6 x 6
#>      V1    V2    V3    V4    V5    V6
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1   -66   -99   -77     7     2   -99
#> 2   -77     2     4     8     8   -99
#> 3     5   -99     7   -66     6     6
#> 4   -99   -99   -77     2   -66   -99
#> 5     9     7   -77   -77     7     1
#> 6     8   -99     2   -77   -99   -99

# Recode -66, -77, and -99 as NA: 
tb[tb == -66 | tb == -77 | tb == -99 ] <- NA

tb
#> # A tibble: 6 x 6
#>      V1    V2    V3    V4    V5    V6
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1    NA    NA    NA     7     2    NA
#> 2    NA     2     4     8     8    NA
#> 3     5    NA     7    NA     6     6
#> 4    NA    NA    NA     2    NA    NA
#> 5     9     7    NA    NA     7     1
#> 6     8    NA     2    NA    NA    NA

sum(is.na(tb))  # => 19 NA values
#> [1] 19

For more sophisticated ways of dealing with (and visualising) NA values, see the R packages mice (Multivariate Imputation by Chained Equations), VIM (Visualization and Imputation of Missing Values) and Amelia II.

Other unusual values: Unexpected values, outliers, etc.

While dealing with missing values is a routine task, finding other unusual values involves some hard thinking, but can also be fun, as it requires the state of mind of a detective who gradually uncovers facts – and aims to ultimately reveal the truth – by questioning a witness or suspect. In principle, we can detect atypical values and outliers (see Exercise 3 of WPA03 for different definitions) by counting values, computing descriptive statistics, or by checking the distributions of raw values. To illustrate this, let’s inspect the data of AHI_CESD further. Since participants in this study were measured repeatedly over a range of several months, two good questions to start with are:

How many participants are there for each occasion?
How many occasions are there for each participant?

The first question concerns a trend (i.e., development over time) and can easily be answered by using dplyr or ggplot:

# How many participants are there for each occasion? 

# Data:
df <- AHI_CESD

# (1) Table of grouped counts:
# Using dplyr pipe: 
id_by_occ <- df %>% 
  group_by(occasion) %>%
  count()
id_by_occ
#> # A tibble: 6 x 2
#> # Groups:   occasion [6]
#>   occasion     n
#>      <int> <int>
#> 1        0   295
#> 2        1   147
#> 3        2   157
#> 4        3   139
#> 5        4   134
#> 6        5   120

# (2) Plot count of IDs by occasion as bar plot:

# (a) Using raw data as input:
ggplot(df, aes(x = occasion)) +
  geom_bar(fill = seeblau) + 
  labs(title = "Number of participants by occasion (from raw data)") +
  theme_bw()


# (b) Using id_by_occ from (1) as input: 
ggplot(id_by_occ, aes(x = occasion)) +
  geom_bar(aes(y = n), stat = "identity", fill = seeblau) + 
  labs(title = "Number of participants by occasion (from table)") +
  theme_bw()

Note the similarities and differences between these two bar plots. Although they look almost the same (except for the label on the y-axis and their title), they use completely different data as inputs. Which of the 2 plot versions would make it easier to add additional information (like error bars, a line indicating the mean counts of participants, or text labels showing the count values for each category)?

# How many occasions are there for each participant?

# Data:
df <- AHI_CESD

# Graphical solution (with ggplot):
ggplot(df) +
  geom_bar(aes(x = id), fill = seeblau) +
  labs(title = "Number of occasions by participant (raw)") +
  theme_bw()


# This does the trick, but looks messy. 
# A cleaner solution uses 2 steps:

# (1) Table of grouped counts: 
occ_by_id <- df %>%
  group_by(id) %>%
  count()
occ_by_id
#> # A tibble: 295 x 2
#> # Groups:   id [295]
#>       id     n
#>    <int> <int>
#>  1     1     2
#>  2     2     6
#>  3     3     2
#>  4     4     4
#>  5     5     5
#>  6     6     3
#>  7     7     1
#>  8     8     6
#>  9     9     1
#> 10    10     1
#> # ... with 285 more rows

# (2) Plot table:

# (a) re-create previous plot (by now from occ_by_id):
ggplot(occ_by_id) +
  geom_bar(aes(x = id, y = n), stat = "identity", fill = seeblau) +
  labs(title = "Number of occasions by participant (2a)") +
  theme_bw()


# (b) reordering id by count n:
ggplot(occ_by_id) +
  geom_bar(aes(x = reorder(id, n), y = n), stat = "identity", fill = seeblau) +
  labs(title = "Number of occasions by participant (2b)") +
  theme_bw()


# (c) more efficient solution:
occ_by_id_2 <- occ_by_id %>% 
  arrange(n)
occ_by_id_2
#> # A tibble: 295 x 2
#> # Groups:   id [295]
#>       id     n
#>    <int> <int>
#>  1     7     1
#>  2     9     1
#>  3    10     1
#>  4    15     1
#>  5    16     1
#>  6    17     1
#>  7    20     1
#>  8    23     1
#>  9    25     1
#> 10    26     1
#> # ... with 285 more rows

## (+) Add rownames (1:n) as a column:
# occ_by_id_2 <- rownames_to_column(occ_by_id_2)  # rowname is character variable!

## (+) Add a variable row that contains index of current row:
occ_by_id_2$row <- 1:nrow(occ_by_id_2)
occ_by_id_2
#> # A tibble: 295 x 3
#> # Groups:   id [295]
#>       id     n   row
#>    <int> <int> <int>
#>  1     7     1     1
#>  2     9     1     2
#>  3    10     1     3
#>  4    15     1     4
#>  5    16     1     5
#>  6    17     1     6
#>  7    20     1     7
#>  8    23     1     8
#>  9    25     1     9
#> 10    26     1    10
#> # ... with 285 more rows

ggplot(occ_by_id_2) +
  geom_bar(aes(x = row, y = n), stat = "identity", fill = seeblau) +
  labs(title = "Number of occasions by participant (2c)") +
  theme_bw()

Our exploration shows that every participant was measured on at least 1 and at most 6 occasions (check range(occ_by_id$n) to verify this). This raises two additional questions:

Was every participant measured on occasion 0 (i.e., the pre-test)?
Was every participant measured only once per occasion?

The first question may seem a bit picky, but do you really know that nobody showed up late (i.e., missed occasion 0) for the study? Actually, we do already know this, since we counted the number of participants and the number of participants per occasion above:

the study sample contains 295 participants (check nrow(posPsy_p_info)), and
the count of participants per occasion showed a value of 295 for occasion 0 in id_by_occ.

For the record, we testify that:

nrow(posPsy_p_info)  # 295 participants
#> [1] 295
id_by_occ$n[1]       # 295 participants counted (as n) in 1st line of id_by_occ
#> [1] 295
nrow(posPsy_p_info) == id_by_occ$n[1]  # TRUE (qed)
#> [1] TRUE

To answer the second question, we can count the number of lines in df per id and occasion and then check whether any unexpected values (different from 1, i.e., n != 1) occur:

# Table (with dplyr):
id_occ <- df %>%
  group_by(id, occasion) %>%
  count()
id_occ
#> # A tibble: 990 x 3
#> # Groups:   id, occasion [990]
#>       id occasion     n
#>    <int>    <int> <int>
#>  1     1        0     1
#>  2     1        1     1
#>  3     2        0     1
#>  4     2        1     1
#>  5     2        2     1
#>  6     2        3     1
#>  7     2        4     1
#>  8     2        5     1
#>  9     3        0     1
#> 10     3        2     1
#> # ... with 980 more rows

summary(id_occ$occasion)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   0.000   2.000   2.028   4.000   5.000
summary(id_occ$n)  # Max is 2 (not 1)!
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   1.000   1.000   1.000   1.002   1.000   2.000

# Do some occasions occur with other counts than 1?
id_occ %>%
  filter(n != 1)
#> # A tibble: 2 x 3
#> # Groups:   id, occasion [2]
#>      id occasion     n
#>   <int>    <int> <int>
#> 1     8        2     2
#> 2    64        4     2

# => Participants with id of 8 and 64: 
#    2 instances of occasion 2 and 4, respectively.

Importantly, 2 participants (8 and 64) are counted twice for an occasion (2 and 4, respectively).

Compare our counts of id_by_occ with Table 1 (p. 4, of Woodworth et al., 2018), which also shows the number of participants who responded on each of the 6 measurement occasions): As the counts in this table correspond to ours, the repeated instances for some measurement occasions (which could indicate data entry errors, but also be due to large variability in the time inteval between measurements) were not reported in the original analysis (i.e., Table 1). This suggests that the 2 participants with repeated occasions were simply counted and measured twice on one occasion and missing from another one.

Another good question to ask is:

How does the number of elapsed.days correspond to the measurement occasions?

This brings us to a more general principle and corresponding practice that involves different ways of viewing distributions of values.

Viewing distributions

Principle 6: Inspect the distributions of variables.

A simple way to get an initial idea about the range of a variable x is to compute summary(x). However, plotting the distribution of individual variables (e.g., with histograms or frequency polygons) provides a more informative overview over the distributions of values present in the data (e.g., whether the assumptions of statistical tests are met) and can provide valuable hints regarding unusual or unexpected values.

Histograms

A histogram provides a cumulative overview of a variable’s values along one dimension (see our examples in visualizing data). Here, we use a histogram of elapsed.days to address the question:

How are the measurement times (elapsed.days) distributed overall (and relative to the stated times of each occasion)?

## Data:
df <- AHI_CESD

summary(df$elapsed.days)  # provides range information:
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00    0.00   14.79   44.31   90.96  223.82

# Vector of official measurement days [based on Table 1 (p. 4)]: 
occ_days <- c(0, 7, 14, 38, 98, 189)
names(occ_days) <- c("0: pre-test", "1: post-test", "2: 1-week", "3: 2-weeks", "4: 3-months", "5: 6-months")
occ_days
#>  0: pre-test 1: post-test    2: 1-week   3: 2-weeks  4: 3-months 
#>            0            7           14           38           98 
#>  5: 6-months 
#>          189

# (a) Histogram: 
ggplot(df, aes(x = elapsed.days)) +
  geom_histogram(fill = seeblau, binwidth = 1) +
  geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) + 
  labs(title = "Distribution of occasions (a: histogram)") +
  theme_bw()

Note: The first 3 occasions are as expected. However, occasions 4 to 6 appear shifted to the left (i.e., were about 7 days earlier than stated).

Alternative ways of viewing distributions

Alternative ways to plot cumulative distributions of values include frequency polygons, density plots, and rug plots.

# (b) frequency polygon:
ggplot(df, aes(x = elapsed.days)) +
  # geom_histogram(fill = seeblau, binwidth = 1) +
  geom_freqpoly(binwidth = 7, color = seeblau) +
  geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) + 
  labs(title = "Distribution of occasions (b: frequency polygon)") +
  theme_bw()


# (c) density plot:
ggplot(df, aes(x = elapsed.days)) +
  # geom_histogram(fill = seeblau, binwidth = 1) +
  geom_density(fill = seeblau) +
  geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) + 
  labs(title = "Distribution of occasions (c: density plot)") +
  theme_bw()


# (d) rug plot:
ggplot(df, aes(x = elapsed.days)) +
  geom_freqpoly(binwidth = 7, color = seeblau) +
  geom_rug(size = 1, color = "black", alpha = 1/4) +
  geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) + 
  labs(title = "Distribution of occasions (d: frequency polygon with rug plot)") +
  theme_bw()

Filtering values

Once we have detected something noteworthy or strange, we may want to mark or exclude some observations (rows) or variables (columns). As we can easily select rows and filter cases (see dplyr::select and dplyr::filter in the last session), it may be good to create and include some dedicated filter variables in our data. Let’s use the participant data posPsy_p_info to illustrate how we can create filter variables and then filter cases:

# Data:
p_info <- posPsy_p_info
dim(p_info)  # 295 x 6
#> [1] 295   6

In previous exercises (see Exercise 5 of WPA02 and Exercise 4 of WPA03), we answered the question:

What is the age range of participants (overall and by intervention)?

by using both ggplot and dplyr pipes. For instance, we can answer questions about the distribution of age values by plotting histograms:

# Age range:
range(p_info$age)
#> [1] 18 83
summary(p_info$age)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   18.00   34.00   44.00   43.76   52.00   83.00

# Histogramm showing the overall distribution of age 
ggplot(p_info) +
  geom_histogram(mapping = aes(age), binwidth = 2, fill = "gold", col = "black") +
  theme_bw() +
  labs(title = "Distribution of age values (overall)")


# Create 4 histogramms showing the distribution of age by intervention:  
ggplot(p_info) +
  geom_histogram(mapping = aes(age), binwidth = 5, fill = seeblau, col = "black") +
  theme_bw() +
  labs(title = "Distribution of age values by intervention (all data)") +
  facet_grid(.~intervention)

For practice purposes, suppose we only wanted to include participants up to an age of 70 years in some analysis. Rather than dropping these participants from the file, we can introduce a filter variable that is TRUE when some criterion (here: age > 70) is satisfied, and otherwise FALSE:

Principle 7: Use filter variables to identify and select sub-sets of observations and variables.

# How many participants are over 70?
sum(p_info$age > 70)  # 6 people with age > 70
#> [1] 6

# Which ones?
which(p_info$age > 70)  #  51  83 114 155 215 244
#> [1]  51  83 114 155 215 244

# Show their details (using dplyr): 
p_info %>%
  filter(age > 70)
#> # A tibble: 6 x 6
#>      id intervention   sex   age  educ income
#>   <int>        <int> <int> <int> <int>  <int>
#> 1    51            2     1    83     2      2
#> 2    83            4     1    71     2      1
#> 3   114            1     1    71     4      1
#> 4   155            3     2    71     4      1
#> 5   215            4     1    75     3      2
#> 6   244            2     1    77     3      2

# Add a corresponding filter variable to df:
p_info <- p_info %>%
  mutate(over_70 = (age > 70))

dim(p_info) # => 7 variables (now including over_70 as a logical variable): 
#> [1] 295   7
head(p_info)
#> # A tibble: 6 x 7
#>      id intervention   sex   age  educ income over_70
#>   <int>        <int> <int> <int> <int>  <int> <lgl>  
#> 1     1            4     2    35     5      3 FALSE  
#> 2     2            1     1    59     1      1 FALSE  
#> 3     3            4     1    51     4      3 FALSE  
#> 4     4            3     1    50     5      2 FALSE  
#> 5     5            2     2    58     5      2 FALSE  
#> 6     6            1     1    31     5      1 FALSE

# Check details again (but applying filter to over_70):
p_info %>%
  filter(over_70)
#> # A tibble: 6 x 7
#>      id intervention   sex   age  educ income over_70
#>   <int>        <int> <int> <int> <int>  <int> <lgl>  
#> 1    51            2     1    83     2      2 TRUE   
#> 2    83            4     1    71     2      1 TRUE   
#> 3   114            1     1    71     4      1 TRUE   
#> 4   155            3     2    71     4      1 TRUE   
#> 5   215            4     1    75     3      2 TRUE   
#> 6   244            2     1    77     3      2 TRUE

In the present case, filter(over_70) is about as long and complicated as filter(age > 70) and thus not really necessary. However, defining explicit filter variables can pay off when constructing more complex filter conditions or when needing several sub-sets of the data (e.g., for cross-validation purposes). Given an explicit filter variable, we can later filter any analysis (or plot) on the fly:

# Age distribution of participants up to 70 (i.e., not over 70):
p_info %>% 
  filter(over_70 == FALSE) %>%
  ggplot() +
  geom_histogram(mapping = aes(age), binwidth = 5, fill = "olivedrab3", col = "black") +
  theme_bw() +
  labs(title = "Distribution of age values by intervention (without participants over 70)") +
  facet_grid(.~intervention)


# Alternatively, we can quickly create sub-sets of the data:
p_info_young <- p_info %>% 
  filter(over_70 == FALSE)

p_info_old <- p_info %>% 
  filter(over_70 == TRUE)

dim(p_info_young)  # 289 participants (295 - 6) x 7 variables 
#> [1] 289   7
dim(p_info_old)    #   6 participants           x 7 variables 
#> [1] 6 7

Viewing relationships

Most research hypotheses involve relationships or measures of correspondence between 2 or more (continuous or categorical) variables. This brings us to another principle:

Principle 8: Inspect relationships between variables.

Scatterplots

A scatterplot shows the relationship between 2 (typically continuous) variables:

# Data:
df <- AHI_CESD
dim(df)  # 992 50
#> [1] 992  50
# df

# Scatterplot (overall):
ggplot(df) +
  geom_point(aes(x = ahiTotal, y = cesdTotal), size = 2, alpha = 1/4) + 
  geom_abline(intercept = 100, slope = -1, col = seeblau) +
  labs(title = "Relationship between ahiTotal and cesdTotal (overall)", 
       x = "ahiTotal", y = "cesdTotal") +
  theme_bw()


# Scatterplot (with facets):
ggplot(df) +
  geom_point(aes(x = ahiTotal, y = cesdTotal), size = 1, alpha = 1/3) + 
  geom_abline(intercept = 100, slope = -1, col = seeblau) +
  labs(title = "Relationship between ahiTotal and cesdTotal (by intervention)", 
       x = "ahiTotal", y = "cesdTotal") +
  facet_wrap(~intervention) + 
  theme_bw()

Bar plots and line graphs

A relationship often of interest in psychology and other sciences are trends or developments over time. When studying trends or developments, the time variable is not always explicitly expressed in units of time (e.g., days, weeks, months), but often encoded implicitly (e.g., as repeated measurements).

Principle 9: Inspect trends over time or repeated measurements.

Trends can be expressed in many ways. Any data that measures some variable more than once and uses some unit of time (e.g., in AHI_CESD we have a choice between occasion vs. elapsed.days) can be used to address the question:

How do values change over time?

To visualize trends, we typically map the temporal variable (in units of time or the counts of repeated measurements) to the x-axis of a plot. If this is the case, many different types of plots can express developments over time. In our session on in visualizing data, we have already encountered bar plots and line graphs, which are well-suited to plot trends. Actually, our bar plots above (using either the raw data of AHI_CESD or our summary table id_by_occ to show the number of participants by occasion) expressed a trend: How many participants dropped out of the study from occasion 0 to 6. We can make the same point in a line plot:

# Data:
df <- AHI_CESD
dim(df)  # 992 50
#> [1] 992  50
df
#> # A tibble: 992 x 50
#>       id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04 ahi05
#>    <int>    <int>        <dbl>        <int> <int> <int> <int> <int> <int>
#>  1     1        0         0               4     2     3     2     3     3
#>  2     1        1        11.8             4     3     3     4     3     3
#>  3     2        0         0               1     3     4     3     4     2
#>  4     2        1         8.02            1     3     4     4     4     3
#>  5     2        2        14.3             1     3     4     4     4     3
#>  6     2        3        32.0             1     3     4     4     4     4
#>  7     2        4        92.2             1     3     3     2     3     3
#>  8     2        5       182.              1     3     3     3     4     2
#>  9     3        0         0               4     3     3     2     4     2
#> 10     3        2        16.4             4     3     3     3     4     4
#> # ... with 982 more rows, and 41 more variables: ahi06 <int>, ahi07 <int>,
#> #   ahi08 <int>, ahi09 <int>, ahi10 <int>, ahi11 <int>, ahi12 <int>,
#> #   ahi13 <int>, ahi14 <int>, ahi15 <int>, ahi16 <int>, ahi17 <int>,
#> #   ahi18 <int>, ahi19 <int>, ahi20 <int>, ahi21 <int>, ahi22 <int>,
#> #   ahi23 <int>, ahi24 <int>, cesd01 <int>, cesd02 <int>, cesd03 <int>,
#> #   cesd04 <int>, cesd05 <int>, cesd06 <int>, cesd07 <int>, cesd08 <int>,
#> #   cesd09 <int>, cesd10 <int>, cesd11 <int>, cesd12 <int>, cesd13 <int>,
#> #   cesd14 <int>, cesd15 <int>, cesd16 <int>, cesd17 <int>, cesd18 <int>,
#> #   cesd19 <int>, cesd20 <int>, ahiTotal <int>, cesdTotal <int>

knitr::kable(id_by_occ)

occasion	n
0	295
1	147
2	157
3	139
4	134
5	120


# Plot count of IDs by occasion as line graph:

# Line plot: Using id_by_occ (from above) as input: 
line_dropouts <- ggplot(id_by_occ, aes(x = occasion)) +
  geom_line(aes(y = n), size = 1.5, color = "gold") + 
  geom_point(aes(y = n), shape = 21, size = 4, stroke = 1, color = "black", fill = "gold") + 
  labs(title = "Number of participants by occasion (line)") +
  theme_bw()
line_dropouts

When comparing this plot with our corresponding bar plot above, we note that ggplot has automatically scaled the y-axis to a range around our frequency values n (here: from about 110 to 300). As a consequence, the number of dropouts over occasions looks more dramatic in the line graph than it did in the bar plot. We can correct for this difference by explicitly scaling the y-axis of the line graph or by plotting both geoms in one plot. And whenever combining multiple geoms in the same plot, we can move common aesthetic elements (e.g., y = n) to the 1st line of our ggplot command, so that we don’t have to repeat them and can only list the specific aesthetics of each geom later:

## Data:
# knitr::kable(id_by_occ)

# Line plot (from above) with corrected y-scale:
line_dropouts + 
  scale_y_continuous(limits = c(0, 300))


# Multiple geoms: 
# Bar + line plot: Using id_by_occ (from above) as input:  
ggplot(id_by_occ, aes(x = occasion, y = n)) +
  geom_bar(stat = "identity", fill = seeblau) + 
  geom_line(size = 1.5, color = "gold") + 
  geom_point(shape = 21, size = 4, stroke = 1, color = "black", fill = "gold") + 
  labs(title = "Number of participants by occasion (bar + line)") +
  theme_bw()

Two key questions in the context of our study are:

How do participants’ scores of happiness (ahiTotal) and depression (cesdTotal) change over time?
Do these trends over time vary by a participant’s type of intervention?

Visualizing these trends will take us more than half-way towards our conclusion, except for a few statistical tests. To answer the first question for the happiness scores (ahiTotal), we compute our desired values with dplyr and then plot the resulting table as a line graph (with an automatic and a full y-scale):

# Data:
df <- AHI_CESD
dim(df)  # 992 50
#> [1] 992  50
# df

# Table 1: 
# Frequency n, mean ahiTotal, and mean cesdTotal by occasion: 
tab_scores_by_occ <- df %>%
  group_by(occasion) %>%
  summarise(n = n(),
            ahiTotal_mn = mean(ahiTotal), 
            cesdTotal_mn = mean(cesdTotal)
            )
knitr::kable(head(tab_scores_by_occ))

occasion	n	ahiTotal_mn	cesdTotal_mn
0	295	69.70508	15.06441
1	147	72.35374	12.96599
2	157	73.04459	12.36943
3	139	74.55396	11.74101
4	134	75.19403	11.70896
5	120	75.83333	12.83333


# Line graph of scores over occasions:
plot_ahi_trend <- ggplot(tab_scores_by_occ, aes(x = occasion, y = ahiTotal_mn)) +
  # geom_bar(aes(y = n), stat = "identity", fill = seeblau) + 
  geom_line(size = 1.5, color = "forestgreen") + 
  geom_point(shape = 21, size = 4, stroke = 1, color = "black", fill = "forestgreen") + 
  labs(title = "ahiTotal_mn by occasion (line)") +
  theme_bw()
plot_ahi_trend


# Line graph with corrected y-scale:
plot_ahi_trend + 
  scale_y_continuous(limits = c(0, 80))

Practice: Create the same plots for mean depression scores cesdTotal. Can you plot both lines in the same plot? What do you see?

To study the trends of scores by participant’s type of intervention, we need to slightly adjust our workflow by adding the intervention variable into our summary table, so that we can later use it in our plotting commands. As before, we will take care of happiness (ahiTotal) and leave it to you to practice on depression (cesdTotal):

# Data:
df <- AHI_CESD
dim(df)  # 992 50
#> [1] 992  50
# df

# Table 2: 
# Scores (n, ahiTotal_mn, and cesdTotal_mn) by occasion AND intervention: 
tab_scores_by_occ_iv <- df %>%
  group_by(occasion, intervention) %>%
  summarise(n = n(),
            ahiTotal_mn = mean(ahiTotal), 
            cesdTotal_mn = mean(cesdTotal)
            )
dim(tab_scores_by_occ_iv)  # 24 x 5
#> [1] 24  5

knitr::kable(head(tab_scores_by_occ_iv))  # print table

occasion	intervention	n	ahiTotal_mn	cesdTotal_mn
0	1	72	68.38889	15.05556
0	2	76	68.75000	16.21053
0	3	74	70.18919	16.08108
0	4	73	71.50685	12.84932
1	1	30	69.53333	15.30000
1	2	48	71.58333	14.58333


# Line graph of scores over occasions AND interventions:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = ahiTotal_mn, 
                                 group = intervention)) +
  geom_line(size = 1.5, color = "forestgreen") + 
  geom_point(shape = 21, size = 4, stroke = 1, color = "black", fill = "forestgreen") + 
  labs(title = "ahiTotal_mn by occasion and intervention (a: lines + points)") +
  theme_bw()

Note that we added group = intervention to the general aesthetics to instruct geom_line to interpret instances of the same intervention but different occasions as belonging together (i.e., to the same line). When using multiple lines, our previous color settings (in geom_line and geom_point) no longer make sense. Instead, we would like to vary the color of our lines and points by intervention. We can achieve this by also moving these aesthetics to the general aes() of the plot:

# Line graph of scores over occasions AND interventions:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = ahiTotal_mn, 
                                 group = intervention, 
                                 color = intervention, 
                                 fill = intervention)) +
  geom_line(size = 1.5) + 
  geom_point(shape = 21, size = 4, stroke = 1, color = "black") + 
  labs(title = "ahiTotal_mn by occasion and intervention (b: lines + points)") +
  theme_bw()

The resulting plot is still suboptimal. The automatic choice of a color scale (in “shades of blue”) and the corresponding legend (on the right) indicates that intervention (i.e., a variable of type integer in df and tab_scores_by_occ_iv) was interpreted as a continuous variable. To turn it into a discrete variable, we can use factor(intervention) for both the color and the fill aesthetics:

# Line graph of scores over occasions AND interventions:
plot_ahi_by_occ_iv <- ggplot(tab_scores_by_occ_iv, 
                             aes(x = occasion, y = ahiTotal_mn, 
                                 group = intervention, 
                                 color = factor(intervention), 
                                 fill = factor(intervention))) +
  geom_line(size = 1.5) + 
  geom_point(shape = 21, size = 4, stroke = 1, color = "black") + 
  labs(title = "ahiTotal_mn by occasion and intervention (c: lines + points, IVs by color)") +
  theme_bw()
plot_ahi_by_occ_iv

An more honest version of the same plot could manually set a different color scheme, but would still look a lot more sobering:

# Line graph of scores over occasions AND interventions:
plot_ahi_by_occ_iv + 
  scale_y_continuous(limits = c(0, 80)) + 
  scale_color_brewer(palette = "Set1") + 
  scale_fill_brewer(palette = "Set1") + 
  labs(title = "Mean ahiTotal scores by occasion and intervention (d)",
       x = "Occasion", y = "Happiness (mean ahiTotal)", 
       color = "Intervention:", fill = "Intervention:")

Practice: Create the same plots for mean depression scores cesdTotal. What do you see?

When viewing and contrasting the trends of groups, we can always zoom in further:

How does this look for individual participants?

Well, how nice of you to ask – so let’s see: We can plot individual trends by occasion or by elapsed days. As there are 295 participants, we use transparency and different facets (by intervention) to help de-cluttering the plot:

# Data:
df <- AHI_CESD
dim(df)  # 992 50
#> [1] 992  50
df
#> # A tibble: 992 x 50
#>       id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04 ahi05
#>    <int>    <int>        <dbl>        <int> <int> <int> <int> <int> <int>
#>  1     1        0         0               4     2     3     2     3     3
#>  2     1        1        11.8             4     3     3     4     3     3
#>  3     2        0         0               1     3     4     3     4     2
#>  4     2        1         8.02            1     3     4     4     4     3
#>  5     2        2        14.3             1     3     4     4     4     3
#>  6     2        3        32.0             1     3     4     4     4     4
#>  7     2        4        92.2             1     3     3     2     3     3
#>  8     2        5       182.              1     3     3     3     4     2
#>  9     3        0         0               4     3     3     2     4     2
#> 10     3        2        16.4             4     3     3     3     4     4
#> # ... with 982 more rows, and 41 more variables: ahi06 <int>, ahi07 <int>,
#> #   ahi08 <int>, ahi09 <int>, ahi10 <int>, ahi11 <int>, ahi12 <int>,
#> #   ahi13 <int>, ahi14 <int>, ahi15 <int>, ahi16 <int>, ahi17 <int>,
#> #   ahi18 <int>, ahi19 <int>, ahi20 <int>, ahi21 <int>, ahi22 <int>,
#> #   ahi23 <int>, ahi24 <int>, cesd01 <int>, cesd02 <int>, cesd03 <int>,
#> #   cesd04 <int>, cesd05 <int>, cesd06 <int>, cesd07 <int>, cesd08 <int>,
#> #   cesd09 <int>, cesd10 <int>, cesd11 <int>, cesd12 <int>, cesd13 <int>,
#> #   cesd14 <int>, cesd15 <int>, cesd16 <int>, cesd17 <int>, cesd18 <int>,
#> #   cesd19 <int>, cesd20 <int>, ahiTotal <int>, cesdTotal <int>

# Individual trends of ahiTotal by occasion and intervention: 
ggplot(df, aes(x = occasion, y = ahiTotal, group = id, color = factor(intervention))) +
  geom_line(size = .5, alpha = .33) + 
  facet_wrap(~intervention) +
  labs(title = "Individual's ahiTotal by occasion for each intervention", 
       color = "Intervention:") +
  scale_color_brewer(palette = "Set1") +
  theme_bw()


# Individual trends of ahiTotal by elapsed.days and intervention: 
ggplot(df, aes(x = elapsed.days, y = ahiTotal, group = id, color = factor(intervention))) +
  geom_line(size = .5, alpha = .33) + 
  facet_wrap(~intervention) +
  labs(title = "Individual's ahiTotal by elapsed.days for each intervention",
       color = "Intervention:") +
  scale_color_brewer(palette = "Set1") +
  theme_bw()

While it is difficult to detect any clear pattern in these plots, it is possible to see some interesting cases that may require closer scrutiny.

Practice: Create the same plots for the individual trends in depression scores cesdTotal. What do you see?

Jitter, box, and violin plots

A measure of correspondence that is common in psychology asks whether the values of some continuous variable vary as a function of the levels of a categorical variable. With regard to our data in AHI_CESD we may wonder:

Do the happiness scores (ahiTotal) vary by intervention?

To visualize the relationship, we cannot use a scatterplot with the mapping x = intervention and y = ahiTotal, as there would only be 4 distinct values for x (go ahead plotting it, if you want to see it). Fortunately, there’s a range of alternatives that allow plotting the raw values and distributions of a continuous variable as a function of a categorical one:

## Data:
# df <- AHI_CESD
# dim(df)  # 992 50
# df

# (a) Jitterplot:
ggplot(df) +
  geom_jitter(aes(x = intervention, y = ahiTotal), width = .1, size = 2, alpha = 1/4) + 
  labs(title = "Values of ahiTotal by intervention (a: jitter)", 
       x = "intervention", y = "ahiTotal") +
  theme_bw()


# (b) Box plot:
ggplot(df) +
  geom_boxplot(aes(x = factor(intervention), y = ahiTotal), color = seeblau, fill = "grey95") + 
  labs(title = "Values of ahiTotal by intervention (b: boxplot)", 
       x = "intervention", y = "ahiTotal") +
  theme_bw()


# Note that we use factor(intervention), as intervention is an integer (i.e., contiuous) variable.

# (c) Violin plot:
ggplot(df) +
  geom_violin(aes(x = factor(intervention), y = ahiTotal), size = 1.5, color = seeblau, fill = "whitesmoke") + 
  labs(title = "Values of ahiTotal by intervention (c: violin)", 
       x = "intervention", y = "ahiTotal") +
  theme_bw()

Note that we sometimes use factor(intervention), as intervention in is an integer (i.e., contiuous) variable.

Sometimes combining 2 or more geoms can be more informative than just using one:

# (d) Combining jitter with boxplot:
ggplot(df) +
  geom_boxplot(aes(x = factor(intervention), y = ahiTotal), color = seeblau, fill = "grey95", alpha = 1) + 
  geom_jitter(aes(x = intervention, y = ahiTotal), width = .1, size = 2, alpha = 1/4) + 
  labs(title = "Values of ahiTotal by intervention (d: jittered boxes)", 
       x = "intervention", y = "ahiTotal") +
  theme_bw()


# (e) Combining violin plot with jitter:
ggplot(df) +
  geom_violin(aes(x = factor(intervention), y = ahiTotal), size = 1.5, color = seeblau) +
  geom_jitter(aes(x = intervention, y = ahiTotal), width = .1, size = 2, alpha = 1/4) + 
  labs(title = "Values of ahiTotal by intervention (e: jittered violins)", 
       x = "intervention", y = "ahiTotal") +
  theme_bw()


# (f) Combining violin plot with boxplot and jitter:
ggplot(df, aes(x = factor(intervention), y = ahiTotal, color = factor(intervention))) +
  geom_violin(size = 1.5, color = seeblau) +
  geom_boxplot(width = .30, color = "grey20", fill = "grey90") + 
  geom_jitter(width = .05, size = 2, alpha = 1/3) + 
  labs(title = "Values of ahiTotal by intervention (f: jittered violins with boxes)", 
       x = "intervention", y = "ahiTotal", color = "Intervention:") + 
  scale_color_brewer(palette = "Set1") +
  theme_bw()

When using multiple geoms in one plot, their order matters, as later geoms are printed on top of earlier ones. In addition, combining geoms typically requires playing with aesthetics. Incidentally, note how the last plot moved some redundant aesthetic mappings (i.e., x, y, and color) from the geoms to the 1st line of the command (i.e., from the mapping argument of the geoms to the general mapping argument).

Tile plots

We already discovered that bar plots can either count cases or show some pre-computed values (when using stat = "identity", see Bar plots in Visualizing data. In the following, we show that we can also compute summary tables and then display their values by the color fill gradient of tiles, or by the size of points:

# Frequency of observations by occasion (x) and intervention (y):

# Data:
df <- AHI_CESD
# dim(df)  # 992 50
# df

# Table 1: 
# Frequency n, mean ahiTotal, and mean cesdTotal by occasion: 
# Use tab_scores_by_occ (from above): 
knitr::kable(head(tab_scores_by_occ))

occasion	n	ahiTotal_mn	cesdTotal_mn
0	295	69.70508	15.06441
1	147	72.35374	12.96599
2	157	73.04459	12.36943
3	139	74.55396	11.74101
4	134	75.19403	11.70896
5	120	75.83333	12.83333


# Table 2: 
# Scores (n, ahiTotal_mn, and cesdTotal_mn) by occasion AND intervention: 
# Use tab_scores_by_occ_iv (from above):
knitr::kable(head(tab_scores_by_occ_iv))

occasion	intervention	n	ahiTotal_mn	cesdTotal_mn
0	1	72	68.38889	15.05556
0	2	76	68.75000	16.21053
0	3	74	70.18919	16.08108
0	4	73	71.50685	12.84932
1	1	30	69.53333	15.30000
1	2	48	71.58333	14.58333


# Tile plot 1: 
# Frequency (n) of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
  geom_tile(aes(fill = n), color = "white", size = .1) + 
  geom_text(aes(label = n), color = "white") +
  labs(title = "Frequency of each intervention by occasion (tile)", 
       x = "Occasion", y = "Intervention", fill = "Frequency:") + 
  theme_bw()


# Tile plot 2: 
# ahiTotal_mn of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
  geom_tile(aes(fill = ahiTotal_mn), color = "grey75", size = 1) + 
  labs(title = "Mean ahiTotal score of each intervention by occasion (tile)", 
       x = "Occasion", y = "Intervention", fill = "ahiTotal_mn:") + 
  scale_fill_gradient(low = "black", high = "gold") +
  theme_bw()

Note the fancy color gradient of the last tile plot: The lowest average scores of happiness (or ahiTotal_mn) are shown in "black", the highest ones in "gold". As the scores vary on a continuous scale, using a continuous color scale is appropriate here. And assigning high happiness to “gold” hopefully makes the plot easier to interpret.

Practice: Create a tile plot that shows the developments of mean depression scores cesdTotal by measurement occasion and intervention. Choose a color scale that you find appropriate. What do you see?

Hint: Your result could look like the following plot:

Point size plots

Another way of expressing the value of a continuous variable – like a frequency n or an average score ahiTotal_mn – as a function of 2 categorical variables is by mapping it to the size dimension of a point. In the following, we illustrate this by creating the same plot twice: First from raw data (df, using geom_count) and then from our summary table (tab_scores_by_occ_iv from above, using geom_point with size mapped to our frequency count n):

# Data:
df <- AHI_CESD
# dim(df)  # 992 50
# df

# Point plot 1a:
# Count frequency of values for each interventin by occasion:
ggplot(df, aes(x = occasion, y = intervention)) +
  geom_count(color = seeblau) + 
  labs(title = "Frequency of each intervention by occasion (count)", 
       x = "Occasion", y = "Intervention", size = "Count:") + 
  scale_size_area(max_size = 12) +
  scale_y_continuous(limits = c(.5, 4.5)) +
  theme_bw()


# Point plot 1b:
# Frequency (n mapped to size) of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
  geom_point(aes(size = n), color = seeblau) + 
  labs(title = "Frequency of each intervention by occasion (point size)", 
       x = "Occasion", y = "Intervention", size = "Frequency:") + 
  scale_size_area(max_size = 12) +
  scale_y_continuous(limits = c(.5, 4.5)) +
  theme_bw()

Tuning plots

As the examples above have shown, most default plots can be modified – and ideally improved – by fine-tuning their visual appearance. In ggplot2, this means setting different aesthetics (or aes() arguments) and scale options (which take different arguments, depending on the type of aesthetic involved). Popular levers for prettifying your plots include:

colors: can be set by the color or fill arguments of geoms (variable when inside aes(...), fixed outside), or by choosing or designing specific color scales (with scale_color);
labels: labs(...) allows setting titles, captions, axis and legend labels, etc.;
legends: can be (re-)moved or edited;
themes: can be selected or modified.

Mixing geoms and aesthetics

By combining the geom_tile and geom_point plots from above with geom_text(aes(label = n) we could simultaneously express 2 different continuous variables on 2 different dimensions (here fill color or size vs. the value shown by a label of text):

# Combo plot 1: 
# ahiTotal_mn of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
  geom_tile(aes(fill = ahiTotal_mn), color = "grey95", size = 1) + 
  geom_text(aes(label = n), color = "white", size = 5, fontface = 2) +
  labs(title = "Mean ahiTotal score of each intervention by occasion (tile)", 
       x = "Occasion", y = "Intervention", fill = "ahiTotal_mn:") + 
  scale_fill_gradient(low = "grey30", high = "gold") +
  theme_bw()


# Combo plot 2b:
# ahiTotal_mn (mapped to point size) and n (text label) for each intervention x occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
  geom_point(aes(size = ahiTotal_mn), color = "forestgreen") + 
  geom_text(aes(label = n), color = "white", size = 3, fontface = 2) +
  labs(title = "ahiTotal_mn and n by intervention x occasion (point size + text label)", 
       x = "Occasion", y = "Intervention", size = "Frequency:") + 
  scale_size_continuous(range = c(5, 11)) +
  scale_y_continuous(limits = c(.5, 4.5)) +
  theme_bw()

However, this mix of multiple dimensions tends to get too confusing, even for ggplot enthusiasts. So let’s conclude by a cautionary note.

Plotting with a purpose

Our last examples have shown that we should better create 2 separate graphics when making multiple points or expressing more than a single relationship in a plot. This brings us to:

Principle 10: A plot should convey its message as clearly as possible.

So go ahead and use the awesome ggplot machine to combine multiple geoms and adjust their aesthetics for as long as this makes your plot prettier. However, rather than getting carried away by the plethora of options, always keep in mind the message that you want to convey. If additional fiddling with aesthetics does not help to clarify your point, you are probably wasting your time by trying to do too much.

10 principles of EDA

Here is a list of the 10 principles covered above:

Start with a clean slate and explicitly load all data and all packages required.
Structure and comment your analysis.
Make copies (and copies of copies) of your data.
Know your data (variables and observations).
Know and deal with unusual values.
Inspect the distributions of variables.
Use filter variables to identify and select sub-sets of observations and variables.
Inspect relationships between variables.
Inspect trends over time or repeated measurements.
A plot should convey its message as clearly as possible.

Adhering to these principles does not guarantee an interesting result, but minimizes the risk of missing something important or violating key assumptions of statistical tests.

Exercises (WPA04)

Working through the above examples should have provided you with a pretty clear picture of the data and results of AHI_CESD. However, we haven’t touched the data contained in posPsy_wide yet. Given that this file was described as a transformed and corrected version of AHI_CESD, we should not expect to find completely different results in it. So let’s use posPsy_wide in some exercises to improve our skills in EDA.

# Load packages:
library(tidyverse)
library(knitr)

# Load data:
# 4. Transformed and corrected version of all data (in wide format): 
posPsy_wide <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_data_wide.csv")

Exercises (WPA04)

Exercise 1: Explore `posPsy_wide`

Load the data stored as http://rpository.com/ds4psy/data/posPsy_data_wide.csv into posPsy_wide (if you haven’t done so above) and explore it.

What are the data’s dimensions, variables, and types of variables? How would you describe the structure of this dataset?

# Load packages:
library(tidyverse)

# Load data:
# 4. Transformed and corrected version of all data (in wide format): 
posPsy_wide <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_data_wide.csv")

df <- posPsy_wide  # copy data

## Checking: 
df
#> # A tibble: 295 x 294
#>       id intervention   sex   age  educ income occasion.0 elapsed.days.0
#>    <int>        <int> <int> <int> <int>  <int>      <int>          <int>
#>  1     1            4     2    35     5      3          0              0
#>  2     2            1     1    59     1      1          0              0
#>  3     3            4     1    51     4      3          0              0
#>  4     4            3     1    50     5      2          0              0
#>  5     5            2     2    58     5      2          0              0
#>  6     6            1     1    31     5      1          0              0
#>  7     7            3     1    44     5      2          0              0
#>  8     8            2     1    57     4      2          0              0
#>  9     9            1     1    36     4      3          0              0
#> 10    10            2     1    45     4      3          0              0
#> # ... with 285 more rows, and 286 more variables: ahi01.0 <int>,
#> #   ahi02.0 <int>, ahi03.0 <int>, ahi04.0 <int>, ahi05.0 <int>,
#> #   ahi06.0 <int>, ahi07.0 <int>, ahi08.0 <int>, ahi09.0 <int>,
#> #   ahi10.0 <int>, ahi11.0 <int>, ahi12.0 <int>, ahi13.0 <int>,
#> #   ahi14.0 <int>, ahi15.0 <int>, ahi16.0 <int>, ahi17.0 <int>,
#> #   ahi18.0 <int>, ahi19.0 <int>, ahi20.0 <int>, ahi21.0 <int>,
#> #   ahi22.0 <int>, ahi23.0 <int>, ahi24.0 <int>, cesd01.0 <int>,
#> #   cesd02.0 <int>, cesd03.0 <int>, cesd04.0 <int>, cesd05.0 <int>,
#> #   cesd06.0 <int>, cesd07.0 <int>, cesd08.0 <int>, cesd09.0 <int>,
#> #   cesd10.0 <int>, cesd11.0 <int>, cesd12.0 <int>, cesd13.0 <int>,
#> #   cesd14.0 <int>, cesd15.0 <int>, cesd16.0 <int>, cesd17.0 <int>,
#> #   cesd18.0 <int>, cesd19.0 <int>, cesd20.0 <int>, ahiTotal.0 <int>,
#> #   cesdTotal.0 <int>, occasion.1 <int>, elapsed.days.1 <dbl>,
#> #   ahi01.1 <int>, ahi02.1 <int>, ahi03.1 <int>, ahi04.1 <int>,
#> #   ahi05.1 <int>, ahi06.1 <int>, ahi07.1 <int>, ahi08.1 <int>,
#> #   ahi09.1 <int>, ahi10.1 <int>, ahi11.1 <int>, ahi12.1 <int>,
#> #   ahi13.1 <int>, ahi14.1 <int>, ahi15.1 <int>, ahi16.1 <int>,
#> #   ahi17.1 <int>, ahi18.1 <int>, ahi19.1 <int>, ahi20.1 <int>,
#> #   ahi21.1 <int>, ahi22.1 <int>, ahi23.1 <int>, ahi24.1 <int>,
#> #   cesd01.1 <int>, cesd02.1 <int>, cesd03.1 <int>, cesd04.1 <int>,
#> #   cesd05.1 <int>, cesd06.1 <int>, cesd07.1 <int>, cesd08.1 <int>,
#> #   cesd09.1 <int>, cesd10.1 <int>, cesd11.1 <int>, cesd12.1 <int>,
#> #   cesd13.1 <int>, cesd14.1 <int>, cesd15.1 <int>, cesd16.1 <int>,
#> #   cesd17.1 <int>, cesd18.1 <int>, cesd19.1 <int>, cesd20.1 <int>,
#> #   ahiTotal.1 <int>, cesdTotal.1 <int>, occasion.2 <int>,
#> #   elapsed.days.2 <dbl>, ahi01.2 <int>, ahi02.2 <int>, ahi03.2 <int>,
#> #   ahi04.2 <int>, …
dim(df)  # 295 observations (participants) x 294 variables
#> [1] 295 294
# names(df)
# glimpse(df)

Structure:

posPsy_wide contains 295 observations (participants) and 294 variables.
The first 6 variables (from id to income) describe each particpant, while the rest contains 6 blocks of 48 variables (from occasion.i to cesdTotal.i, and $i = 0, 1,... 5$). [Check: $6$ variables $+ 6 \cdot\ 48$ variables $= 294$ variables.]
Most variables are numbers of type integer, except for elapsed.days.1 to elapsed.days.5, which are numeric of type double.

Thus, we suspect that the data combines the information contained in posPsy_p_info and AHI_CESD, but – in contrast to AHI_CESD – presents it in wide format (1 row per participant) rather than long format (1 row per measurement).

Verify that posPsy_wide contains the same participants as posPsy_p_info.

# Select the first 6 variables of posPsy_wide:
p_info_2 <- posPsy_wide %>%
  select(id:income)

# Verify that they are equal to posPsy_p_info: 
all.equal(p_info_2, posPsy_p_info)
#> [1] TRUE

Inspect posPsy_wide for missing values. Do you see some systematic pattern to them?

df <- posPsy_wide  # copy data
df
#> # A tibble: 295 x 294
#>       id intervention   sex   age  educ income occasion.0 elapsed.days.0
#>    <int>        <int> <int> <int> <int>  <int>      <int>          <int>
#>  1     1            4     2    35     5      3          0              0
#>  2     2            1     1    59     1      1          0              0
#>  3     3            4     1    51     4      3          0              0
#>  4     4            3     1    50     5      2          0              0
#>  5     5            2     2    58     5      2          0              0
#>  6     6            1     1    31     5      1          0              0
#>  7     7            3     1    44     5      2          0              0
#>  8     8            2     1    57     4      2          0              0
#>  9     9            1     1    36     4      3          0              0
#> 10    10            2     1    45     4      3          0              0
#> # ... with 285 more rows, and 286 more variables: ahi01.0 <int>,
#> #   ahi02.0 <int>, ahi03.0 <int>, ahi04.0 <int>, ahi05.0 <int>,
#> #   ahi06.0 <int>, ahi07.0 <int>, ahi08.0 <int>, ahi09.0 <int>,
#> #   ahi10.0 <int>, ahi11.0 <int>, ahi12.0 <int>, ahi13.0 <int>,
#> #   ahi14.0 <int>, ahi15.0 <int>, ahi16.0 <int>, ahi17.0 <int>,
#> #   ahi18.0 <int>, ahi19.0 <int>, ahi20.0 <int>, ahi21.0 <int>,
#> #   ahi22.0 <int>, ahi23.0 <int>, ahi24.0 <int>, cesd01.0 <int>,
#> #   cesd02.0 <int>, cesd03.0 <int>, cesd04.0 <int>, cesd05.0 <int>,
#> #   cesd06.0 <int>, cesd07.0 <int>, cesd08.0 <int>, cesd09.0 <int>,
#> #   cesd10.0 <int>, cesd11.0 <int>, cesd12.0 <int>, cesd13.0 <int>,
#> #   cesd14.0 <int>, cesd15.0 <int>, cesd16.0 <int>, cesd17.0 <int>,
#> #   cesd18.0 <int>, cesd19.0 <int>, cesd20.0 <int>, ahiTotal.0 <int>,
#> #   cesdTotal.0 <int>, occasion.1 <int>, elapsed.days.1 <dbl>,
#> #   ahi01.1 <int>, ahi02.1 <int>, ahi03.1 <int>, ahi04.1 <int>,
#> #   ahi05.1 <int>, ahi06.1 <int>, ahi07.1 <int>, ahi08.1 <int>,
#> #   ahi09.1 <int>, ahi10.1 <int>, ahi11.1 <int>, ahi12.1 <int>,
#> #   ahi13.1 <int>, ahi14.1 <int>, ahi15.1 <int>, ahi16.1 <int>,
#> #   ahi17.1 <int>, ahi18.1 <int>, ahi19.1 <int>, ahi20.1 <int>,
#> #   ahi21.1 <int>, ahi22.1 <int>, ahi23.1 <int>, ahi24.1 <int>,
#> #   cesd01.1 <int>, cesd02.1 <int>, cesd03.1 <int>, cesd04.1 <int>,
#> #   cesd05.1 <int>, cesd06.1 <int>, cesd07.1 <int>, cesd08.1 <int>,
#> #   cesd09.1 <int>, cesd10.1 <int>, cesd11.1 <int>, cesd12.1 <int>,
#> #   cesd13.1 <int>, cesd14.1 <int>, cesd15.1 <int>, cesd16.1 <int>,
#> #   cesd17.1 <int>, cesd18.1 <int>, cesd19.1 <int>, cesd20.1 <int>,
#> #   ahiTotal.1 <int>, cesdTotal.1 <int>, occasion.2 <int>,
#> #   elapsed.days.2 <dbl>, ahi01.2 <int>, ahi02.2 <int>, ahi03.2 <int>,
#> #   ahi04.2 <int>, …
sum(is.na(df))  # 37440
#> [1] 37440

# Hypothesis 1: ------ 
# When an occasion is missing, all corresponding variables are missing as well.

# Test for occasion.1: 
df %>%
  filter(is.na(occasion.1)) %>%       # occasion.1 is missing: 
  select(occasion.1:cesdTotal.1) %>%  # block of 48 variables corresponding to this occasion
  sum(., na.rm = TRUE)                # => 0 (i.e., only NA values left).
#> [1] 0

# Test for occasion.2:
df %>%
  filter(is.na(occasion.2)) %>%       # occasion.2 is missing: 
  select(occasion.2:cesdTotal.2) %>%  # block of 48 variables corresponding to this occasion
  sum(., na.rm = TRUE)                # => 0 (i.e., only NA values left).
#> [1] 0

# Test for occasion.3:
df %>%
  filter(is.na(occasion.3)) %>%       # occasion.3 is missing: 
  select(occasion.3:cesdTotal.3) %>%  # block of 48 variables corresponding to this occasion
  sum(., na.rm = TRUE)                # => 0 (i.e., only NA values left).
#> [1] 0

# Test for occasion.4:
df %>%
  filter(is.na(occasion.4)) %>%       # occasion.4 is missing: 
  select(occasion.4:cesdTotal.4) %>%  # block of 48 variables corresponding to this occasion
  sum(., na.rm = TRUE)                # => 0 (i.e., only NA values left).
#> [1] 0

# Test for occasion.5:
df %>%
  filter(is.na(occasion.5)) %>%       # occasion.5 is missing: 
  select(occasion.5:cesdTotal.5) %>%  # block of 48 variables corresponding to this occasion
  sum(., na.rm = TRUE)                # => 0 (i.e., only NA values left).
#> [1] 0

# Hypothesis 2: ------ 
# If missing measurements are the ONLY reason for NA values, 
# there should be no other missing values (e.g., individual NA values within measurements).

# Test whether missing occasions account for all NA values:
all_missings <- sum(is.na(df))  # 37440
occ_missings <- sum(is.na(df$occasion.1), is.na(df$occasion.2), is.na(df$occasion.3), 
                    is.na(df$occasion.4), is.na(df$occasion.5)) * 48 
occ_missings == all_missings    # TRUE (qed).
#> [1] TRUE

Exercise 2: Income

Transform the data in posPsy_wide to compute the distribution of income levels and compare your results to those of WPA03: Exercise 4 and WPA01: Exercise 6, which both used the p_info data file http://rpository.com/ds4psy/data/posPsy_participants.csv.

df <- posPsy_wide

# 1. Table showing a distribution of income levels: ------ 
df %>%
  group_by(income) %>%
  summarise(n = n(),
            pc = n/nrow(df) * 100
            )
#> # A tibble: 3 x 3
#>   income     n    pc
#>    <int> <int> <dbl>
#> 1      1    73  24.7
#> 2      2   136  46.1
#> 3      3    86  29.2

# Compare with results from other data set:
p_info <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_participants.csv")

# Using dplyr on p_info data [WPA03 Exercise 4]:
p_info %>%
  group_by(income) %>%
  summarise(n = n(),
            pc = n/nrow(p_info) * 100) 
#> # A tibble: 3 x 3
#>   income     n    pc
#>    <int> <int> <dbl>
#> 1      1    73  24.7
#> 2      2   136  46.1
#> 3      3    86  29.2

# Using p_info and base R [WPA01, Exercise 6]:
n_income_low <- sum(p_info$income < 2)
pc_income_low <- n_income_low/nrow(p_info) * 100
pc_income_low  # 24.74576 
#> [1] 24.74576

# => Results seem to be identical, which is good.

Visualize the distributions of income levels overall and separately for each intervention (Hints: Use dplyr to compute a corresponding summary table before using calling ggplot. To create separate plots, use facet (e.g., facet_wrap(~intervention)) on the overall plot.)

# 2. Visualizations: ------ 

# Plot of income levels overall: 
df %>% 
  group_by(income) %>%
  count() %>%
  ggplot(., aes(fill = factor(income))) + 
    geom_bar(aes(x = income, y = n), stat = "identity") +
    labs(title = "Income levels (a: overall)", 
         x = "Income level", y = "Number", fill = "Income:") + 
    scale_fill_brewer(palette = "Set1") +
    theme_bw()


# Plot of income levels by intervention: 
df %>%
  group_by(intervention, income) %>%
  count() %>%
  ggplot(., aes(fill = factor(income))) + 
    geom_bar(aes(x = income, y = n), stat = "identity") +
    labs(title = "Income levels (b: by intervention)", 
         x = "Income level", y = "Number", fill = "Income:") + 
    scale_fill_brewer(palette = "Set1") +
    theme_bw() +
    facet_wrap(~intervention)

Plot some initial variables (happiness or depression) by income.

+++ here now +++

Exercise +: Total scores?

Use data transformation on the file AHI_CESD to answer the following question:

Are the variables ahiTotal and cesdTotal the sums of all values in ahi01 to ahi24 and cesd01 to cesd20, respectively?

## Data:
## 2. Original DVs in long format:
# AHI_CESD <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_AHI_CESD.csv")

df <- AHI_CESD
dim(df)  # 992 x 50
#> [1] 992  50

# Are the variables `ahiTotal` and `cesdTotal` the sums of all values in `ahi01` to `ahi24` and `cesd01` to `cesd20`, respectively?
# (a) ahi: 
ahi <- df %>% 
  select(ahi01:ahi24)

ahi_sums <- rowSums(ahi)
sum(ahi_sums == df$ahiTotal)  # 992 (qed)
#> [1] 992
all.equal(ahi_sums, df$ahiTotal)  # TRUE
#> [1] TRUE

# (b) cesd:
cesd <- df %>% 
  select(cesd01:cesd20)

cesd_sums <- rowSums(cesd)
cesd_sums
#>  [1] 36 36 27 30 27 26 37 40 35 37 45 51 51 47 43 42 38 29 28 45 54 45 32
#> [24] 46 57 63 59 50 61 45 43 37 35 41 34 36 42 31 45 30 32 50 58 40 40 34
#> [47] 44 42 41 42 42 46 20 36 36 38 34 39 46 43 45 39 41 55 53 42 34 34 37
#> [70] 41 38 47 36 45 31
#>  [ reached getOption("max.print") -- omitted 917 entries ]
sum(cesd_sums == df$cesdTotal)  # 0 => all different!
#> [1] 0

Exercise +: Same plot twice

Plot same plot twice: 1. from a lot of raw data and 2. from much leaner table of aggregated data.

Conclusion

All ds4psy essentials:

Nr.	Topic
0.	Syllabus
1.	Basic R concepts and commands
2.	Visualizing data
3.	Transforming data
4.	Exploring data (EDA)
+.	Datasets

[Last update on 2018-11-29 18:46:50 by hn.]

The participant part of this data has been used as p_info in the previous sessions.↩
Actually, listing everything required to execute a file can get quite excessive (check out sessionInfo() or the session_info and package_info functions of devtools, in case you’re interested, and consider using the packrat package in case you want to preserve or share your current setup). Hence, most people only explicitly list and load non-standard packages (in their current version).↩
RStudio helps with writing well-structured and clean code by providing foldable sections (by typing repeated characters, like ----) or automatic shortcuts (e.g., see the keyboard shortcuts for automatic indenting and un-/commenting).↩

Exploring data (EDA, ds4psy)

Hansjörg Neth, SPDS, uni.kn

2018 11 28

Introduction

Course coordinates

Preparations

Exploring data

Introduction

Key concepts

What is EDA?

Principles and practices

Setting the stage

Cleaning up

Loading packages and data

Clarity

Safety

Copies of objects

Copies of data

Screening data

Basic properties

Unusual values

Missing values

Other unusual values: Unexpected values, outliers, etc.

Viewing distributions

Histograms

Alternative ways of viewing distributions

Filtering values

Viewing relationships

Scatterplots

Bar plots and line graphs

Jitter, box, and violin plots

Tile plots

Point size plots

Tuning plots

Mixing geoms and aesthetics

Plotting with a purpose

10 principles of EDA

Exercises (WPA04)

Exercises (WPA04)

Exercise 1: Explore posPsy_wide

Exercise 2: Income

Exercise +: Total scores?

Exercise +: Same plot twice

More on EDA

Conclusion

Exercise 1: Explore `posPsy_wide`